%% Assignment 10
% For CIE4440

close all;clear all;clc;

%%   Load data
load data

%%  10.1
% Distance matrix

for i=1:length(data)
    xcoor(i)=data{i}.xCoor;
    ycoor(i)=data{i}.yCoor;
end
for i = 1:length(data)
    data{i}.dx=xcoor-xcoor(i);
    data{i}.dy=ycoor-ycoor(i);
    data{i}.distance=sqrt(data{i}.dx.^2+data{i}.dy.^2);
end

% Nominator
for i = 1:length(data)
    data{i}.nominator=NaN(1,length(data));
    for ii=1:length(data)
       data{i}.nominator(ii)= (data{ii}.rainJan-data{i}.rainJan).^2;
    end
end

% sort by distance
distArray=NaN(length(data),length(data));
for i = 1:length(data)
    distArray(i,1:length(data))=data{i}.distance;
    nominatorArray(i,1:length(data))=data{i}.nominator;
    [fromArray{i,1:length(data)}]=deal(repmat(data{i}.name,1,1));
    [toArray{i,1:length(data)}]=deal(data{i}.distancetoname{:});
end
% Only use half the matrix
distArray2=distArray.*triu(ones(length(data),length(data)));

% Sort the matrix
[S,I]=sort(distArray2(:),'ascend');

% Filter the zeros
I=I(S>0);
S=S(S>0);


for i=1:length(I)
    % Find the corresponding nominator and city names
    nomArray(i,1)=nominatorArray(I(i));
    fromArray2{i,1}=fromArray{I(i)};
    toArray2{i,1}=toArray{I(i)};
    
    % Define distance classes
    if S(i)>=0 && S(i)<50, distClass(i)=1; end
    if S(i)>=50 && S(i)<100, distClass(i)=2; end
    if S(i)>=100 && S(i)<150, distClass(i)=3; end
    if S(i)>=150 && S(i)<200, distClass(i)=4; end
    if S(i)>=200 && S(i)<250, distClass(i)=5; end
    if S(i)>=250 && S(i)<300, distClass(i)=6; end
    if S(i)<0 && S(i)>=300, error('define an extra class'); end
end
Array=[num2cell([I,S,distClass',nomArray]),fromArray2,toArray2];

for i=1:6 % six classes defined
    classSize(i)=length([Array{[Array{:,3}]==i,4}]);
    sumPerClass(i)=sum([Array{[Array{:,3}]==i,4}]);
    hAve(i)=i*50-25;
    gamma(i)=.5*sumPerClass(i)/classSize(i);
end

% create monomial line/fittiing curve/variogram model
% omega and lambda are given for a monomial fittine curve
omega=1.17; lambda=2.19;
x=1:300;
y=omega*abs(x).^lambda;

% Plot the calculated values
figure(1)
cax = gca;
set(cax,'NextPlot','add');
scatter(hAve,gamma); 
plot(x,y);
xlabel('h [km]');
ylabel('gamma(h)');

% STEP 7: generate matrix A
% Using Delft instead of Hoogeveen
gammaCalculated=omega.*abs(distArray).^lambda;

% Define Delft location
Delft.name='Delft';
Delft.distancetoname={'Amsterdam','De Bilt','Eindhoven','Den Helder','Groningen','Leeuwarden','Maastricht','Rotterdam','Twente','Vlissingen','Deelen'};
Delft.xCoor=85;
Delft.yCoor=447;
Delft.rainJan=NaN;
Delft.rainJul=NaN;

for i = 1:length(data)
    Delft.dx=xcoor-Delft.xCoor;
    Delft.dy=ycoor-Delft.yCoor;
end
Delft.distance=sqrt(Delft.dx.^2+Delft.dy.^2);

matrixC=omega.*abs(Delft.distance).^lambda;
weight=matrixC/(gammaCalculated);

for i=1:length(data)
   Delft.rainJanMonPart(i)=data{i}.rainJan*weight(i);    
end

% calculate for Delft using Monomial
Delft.rainJan=sum(Delft.rainJanMonPart);
disp(Delft.rainJan);


%% 10.2 find value for Delft
% Already done in 10.1